function thresh = findThreshold(I,options)
    
    acut = 10; %minimum area cutoff - speckles
    aBoundL = 20; %average area in a succesful threshold
    aBoundH = 2000;
    maxa = @(t,I) size(I,1)^2/2-bwarea(I>t);
    t_maxa = fzero(@(t) maxa(t,I),(2^16-1)/2); %get the threshold that gives 1/2 fill
        
    %if c==1
    if 1
        %BF - maximize average area of components
        thresh = fminbnd(@(t) ava(t,I,acut),t_maxa + 3,(2^16-1),options); %padded
        if -ava(thresh,I,acut)<aBoundL || -ava(thresh,I,acut)>aBoundH 
            thresh = fminsearch(@(t) euler(t,I,acut),t_maxa,options)  +  5; %padded to reduce speckling
        end
    else
        %Nuclei - maximize number of components
        thresh = fminbnd(@(t) maxn(t,I,acut),t_maxa,(2^16-1),options);
    end
end

function out = ava(t,I,acut)
cc = bwconncomp(bwareaopen(I>t,acut));
out = -bwarea(bwareaopen(I>t,acut))/cc.NumObjects;
end

function out = areaNum(t,I,acut)
cc = bwconncomp(bwareaopen(I<t,acut));
out = -bwarea(bwareaopen(I<t,acut))*cc.NumObjects^2;
end

function out = maxn(t,I,acut)
cc = bwconncomp(bwareaopen(I>t,acut));
out = -cc.NumObjects; %max number of components
end

function out = euler(t,I,acut)
out = -bweuler(bwareaopen(I>t,acut));
end
